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We solve the impurity problem which arises within nonequilibrium dynamical mean-field theory 
for the Hubbard model by means of a self-consistent perturbation expansion around the atomic limit. 
While the lowest order, known as the non-crossing approximation (NCA), is reliable only when the 
interaction U is much larger than the bandwidth, low-order corrections to the NCA turn out to be 
sufficient to reproduce numerically exact Monte Carlo results in a wide parameter range that covers 
the insulating phase and the metal-insulator crossover regime at not too low temperatures. As an 
application of the perturbative strong-coupling impurity solver we investigate the response of the 
double occupancy in the Mott insulating phase of the Hubbard model to a dynamical change of the 
interaction or the hopping, a technique which has been used as a probe of the Mott insulating state 
in ultracold fermionic gases. 
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I. INTRODUCTION 

Experiments with ultracold atoms in optical lattices,— 
as well as pump-probe spectroscopy with femtosec- 
ond time-resolutioii^ri^ and transport measurements on 
quantum dots^ enable a systematic investigation of 
strongly interacting quantum many-particle systems un- 
der nonequilibrium conditions. In a pump-probe experi- 
ment, a material is excited with a strong laser pulse, and 
its subsequent evolution is probed by a second pulse that 
reaches the sample at a controlled time-delay. The break- 
down of the Mott insulating phase within a few times 
of the inverse hopping has been observed in this way.^^ 
Ultracold gases in optical lattices, on the other hand, 
can be prepared in an equilibrium state and suddenly 
quenched out of equilibrium by modifying a Hamiltonian 
parameter 1^ These experiments can address fundamental 
questions such as the thermalization in isolated quantum 
systems li 

The ongoing experimental progress has stimulated in- 
tensive research on the theoretical side. While many the- 
oretical approaches that are designed for the investiga- 
tion of correlated systems in thermal equilibrium must 
be modified considerably before they can be used to 
compute the real-time evolution, dynamical mean-field 
theory (DMFT)^ is an approximate scheme which is 
per se applicable to both equilibrium and nonequilib- 
rium situations.— The method rehes on a mapping of 
lattice models to a single impurity model, which is ex- 
act in the limit of infinite dimensions ^^ and provides 
a good basis for the realistic simulation of many corre- 
lated materials.— Nonequilibrium DMFT has so far been 
used, e.g., to study transport beyond linear response in 
the Falicov-Kimball model^"— as well as interaction 
quenches and interaction ramps in the Falicov-Kimball 
mode l^^'^^ and in the Hubbard model.— i^ In the last 
part of this paper we will use DMFT to study the re- 
sponse of the Mott insulating phase to a periodic mod- 
ulation of the hopping or the interaction,— similar to 
what can now be done in experiments with cold atomic 



gases i2ii2^ 

Currently, the biggest challenge within the context of 
nonequilibrium DMFT is the development of impurity 
solvers which allow to compute the long-time dynamics 
after a perturbation. An exact solution, via a closed set 
of equations of motion, is known only for the Falicov- 
Kimball model J^i2^ For the Hubbard model, continuous- 
time Quantum Monte Carlo (CTQMC) can in princi- 
ple be used to obtain an unbiased solution.— li^ Both 
the weak-coupling expansior^^i and the strong-couphng 
expansioii^^ of the relevant Anderson impurity model 
have been translated from their respective imaginary- 
time variants (Refs. [2^ and [23 for weak-coupling and 
Ref. |28| for strong-coupling) to the Keldysh formalism, in 
order to study the real-time evolution. However, in these 
real-time Monte Carlo calculations, the accessible times 
are limited by the notorious dynamical sign problem. A 
big advantage of the weak-coupling expansion over the 
strong-coupling expansion is that the diagrammatic se- 
ries simplifies in the case of particle-hole symmetry.— 
On the other hand, the sign problem in a weak-coupling 
calculation increases with the interaction strength. It 
is thus essentially impossible to use CTQMC to study 
complex excitation processes within the Mott insulating 
phase, e.g., the excitation with a short laser pulse and 
the subsequent relaxation. 

In order to avoid the sign problem and access the 
regime of strong interactions and relatively long times, 
we explore the direct summation of the self-consistent 
diagrammatic hybridization expansion up to fixed order, 
as opposed to CTQMC, which is in essence a stochastic 
summation of the full (non self-consistent) series. This 
approach proves to be very accurate in a wide paremeter 
regime and suitable for the calculation of the real-time 
dynamics. 

Systematic approximations for the expansion around 
the atomic limit of the Anderson impurity model have 
been used for a long time.—"— The simplest conserv- 
ing approximation, which has been termed non-crossing 
approximation (NCA), can correctly recover the Kondo 



temperature Tk when charge fluctuations are suppressed 
by the Coulomb interaction [/, although the Fermi-liquid 
behavior for T ^ Tk is not correctly reproduced.— i^ 
If U is finite, however, the width of the Kondo reso- 
nance is severely underestimated, and various resumma- 
tion schemes of the expansion have been devised to cure 
this problem i^^i^^ Already the simplest to NCA within 
these schemes, the so called one-crossing approximation 
(OCA), can cure the deficiencies of NCA to a large ex- 
tent. Motivated by the fact that NCA is already very 
good in the insulating parameter regime. Gull et al^ re- 
cently developed a bold-line hybridization expansion, i.e., 
an approach which is based on a Monte Carlo sampling 
of the corrections beyond NCA. 

Starting with the work of Pruschke, Cox, and Jarrell;^ 
both the NCA and the OCA have been used as an im- 
purity solver for DMFT (for some recent references that 
involve the investigation of real materials, see Refs. WA - 
|44| ). Furthermore, NCA and its corrections can readily be 
translated into the Keldysh formalism to study nonequi- 
librium situations, although the evaluation of higher- 
order diagrams in real time involves quite some numer- 
ical effort. For example, the buildup of the Kondo res- 
onance after a sudden shift of the impurity level in the 
Kondo model has been investigated with NCA4^ The 
fairly accurate results in equilibrium calculations and the 
straightforward portability to the Keldysh contour make 
the self-consistent hybridization expansion an interesting 
candidate for the solution of the impurity problem within 
nonequilibrium DMFT. 

The purpose of this paper is twofold. First, we give 
a detailed description of the self-consistent expansion on 
the Keldysh contour (Sections. Ulland lllll . and we bench- 
mark the method by applying it to the interaction quench 
in the Hubbard model on the Bethe-lattice fSec.lIVp. We 
find similar trends for both equilibrium and nonequilib- 
rium: While NCA is unreliable unless the interaction U 
is much larger than the bandwidth, OCA provides an 
important correction, and the third order is in almost 
quantitative agreement with QMC results over a wide 
parameter range which includes the insulating phase and 
the crossover regime between metal and insulator. As 
an application of the perturbative impurity solver we 
then study the excitation of a Mott insulator by a time- 
dependent modulation of the hopping or the interaction 
strength (Sec.|V|. Because the interaction is rather large 
in this problem, a solution using weak-coupling CTQMC 
is currently not feasible. 



II. DEFINITION OF THE IMPURITY 
PROBLEM ON THE KELDYSH CONTOUR 

To describe a nonequilibrium situation in which the 
system is prepared in a thermal equilibrium state at tem- 
perature T = l/j3 for times t < and later acted on by 
some perturbation, we use the Keldysh formalism4^i^ 
The imaginary-time contour of the Matsubara Green's 
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FIG. 1: The L-shaped contour C for the description of tran- 
sient nonequilibrium states with initial state density matrix 
oc exp[~'l3H{0)]. The indicated time arguments are in cyclic 
order, ii -< f2 ^ is (see text). 



functions for finite-temperature equilibrium states is 
thereby extended to the L-shaped contour C that runs 
from to time imax (i-e., the largest time of interest) 
along the real axis, back to 0, and finally to —i[3 along 
the imaginary time axis (Fig. [1]). The Keldysh formal- 
ism is based on the use of contour-ordered correlation 
functions {TcA{ti)B{t2)), where Tc exchanges the or- 
der of the two operators A{ti) and B{t2) in the product 
A{ti)B{t2) if ^2 appears later on the contour than ii, ac- 
cording to the order which is indicated by the arrows in 
Fig. [TJ An additional minus sign appears if the exchange 
involves an odd number of Fermi operators. The use of 
contour-ordered Green's functions allows the application 
of Wick's theorem if the action is quadratic4^ Depending 
on the choice of the time arguments, a contour-ordered 
correlation function describes either real-time correla- 
tions, or it recovers the imaginary-time ordered corre- 
lation function of the initial equilibrium state. 

In the following sections we consider an impurity model 
which is defined by the following action on the L-shaped 
contour C, 



<5ioc = -i dtHioc[dl{t),dp{t),t], 



(la) 
(lb) 



^hyb = -i dtidt2 y^ dl^jti) Kp^^p^{ti,t2) dpj^t2). (Ic) 

P1-P2 

In this action, dp and dj, denote annihilation and creation 
operators for an electron in the impurity level p {p labels 
spin and orbital degrees of freedom) , and -ffioc is the local 
Hamiltonian of the impurity site, which can be interact- 
ing and time-dependent in general. The hybridization 
function Ap-^^p^(ti,t2) gives the amplitude for the hop- 
ping of an electron from the p2-orbital into the bath at 
time f 2 , its propagation within the bath, and the hopping 
back into the impurity orbital pi at time ii. The action 
P]) can be derived from an impurity Hamiltonian with 
time-dependent coupling between bath and impurity. 



H{t) = Hi,,{t)- 



' / , ^vCjjCy- 



-^[yp,,(i)dtc,-H/i.c.] (2) 



by tracing out the bath degrees of freedom c^. It also 
arises as the effective single-site problem in nonequilib- 
rium DMFT,^^ without direct reference to a given Hamil- 
tonian formulation. 



The single-particle Green's function of the impurity 
model ([l} is given by 



Gp,p'{t,t')^~i{dp{t)dl,{t'))s^,^,, 



(3) 



where the contour-ordered expectation value for the ac- 
tion S is defined as 



{■■■)s = 



Tr[Tc exp(5) ■ ■ ■ 
Tr[Tcexp(5)] 



(4) 



The noncquilibrium formalism presented below reduces 
to the Matsubara formalism for the initial equilibrium 
state when all calculations are restricted to the imaginary 
branch of the contour. For time-arguments t = —it and 
t' = —it' on the imaginary branch of C, the Green's func- 
tion G{t,t') is directly related to the Matsubara Green's 
function C^'It) of the initial thermal equilibrium state, 



G{-iT,~iT') =iG'''{T^T'), 



(5) 



which is translationally invariant in time. An analogous 
equation is used to define the Matsubara component of all 
two-time contour-ordered correlation functions that are 
used in the following. The factor i on the right-hand side 
of Eq. (O is needed to recover the conventional definition 
of the Matsubara functions. 



The local part Hioc{t) of the impurity Hamiltonian is 
diagonalized at each instant of time, 



Hiocit) = J2 \mit))Era{t){m{t)l 



(6) 



and for each eigenstate m one flavor of pseudoparticles, 
with annihilation (creation) operator am , is introduced. 
We assume that eigenstates \m{t)) are smooth functions 
of time. Pseudoparticles are bosons if the state \m) cor- 
responds to an even number of particles on the impu- 
rity, and fermions otherwise. By means of the isomorphy 
\'m{t)) -!->■ ajjjjwac), the physical Hilbert space of the im- 
purity can be identified with the subspace of the pseu- 
doparticle Fock space in which the total number of pseu- 
doparticles, 



/ y ^m'^rm 



(7) 



is exactly one. Hence the expectation value {A{t))s^^ 
of any impurity observable A can be computed in the 
pseudoparticlc space as 

{SqiAU))^ 

(^(0>5,„,p - (^(t))Q=i ^ ,, , ^'°^ (8) 



III. SELF-CONSISTENT DIAGRAMMATIC 

HYBRIDIZATION EXPANSION OF THE 

ANDERSON IMPURITY MODEL 

A. Pseudoparticle representation 

In this section we compute the Green's function ([3]) 
by expanding the expectation value ^ in terms of 
the hybridization function A{t,t'). The self-consistent 
hybridization expansion for the Anderson model in 
thermal equilibrium has been described previously in 
many places r^"—'^'^'^ and the generalization from the 
imaginary-time contour to the Keldysh contour is rather 
straightforward. Nevertheless we give a detailed deriva- 
tion below, in order to discuss some important technical 
differences between the equilibrium and the noncquilib- 
rium variants of the expansion. 

Because the local part (|lb[l of the action is generally 
not quadratic, standard diagrammatic perturbation the- 
ory does not apply to the expansion around the atomic 
limit (A = 0). There exist several related strategies to 
by-pass this difficulty. In the CTQMC variant of the 
hybridization cxpansionj^S high-order time-ordered cor- 
relation functions of the impurity problem are explicitly 
evaluated in a suitable basis of the local problem. We 
will follow a different approach, which is based on the in- 
troduction of auxiliary particlesi^i^ This allows to use 
standard resummation tricks from diagrammatic pertur- 
bation theory, at the expense of having to do a projec- 
tion from the extended pseudoparticle Hilbert space to 
the physical Hilbert space. 



provided that the pseudoparticle action Simp and the ob- 
servable A are constructed such that they coincide with 
Simp and A in the Q = 1 subspace {Sq^i is the projection 
onto Q — 1). The requirement is satisfied by choosing 

4(t)=^F^„(t)aLa„ (9a) 

dp(t)=^F„PJO*aLan (9b) 

F„Ut) = {m{t)\dl\nit)) (9c) 

for the electron annihilation and creation operators, and 



^loc = ^*X] / dt Emit) a\^{t)am{t) 



(10a) 



^hyb — 



E 



m,n.,m' .n' p,p' 



dtdt alj^{t)an{t) x 



X FUt)Kyit,t')Fl'm'it)al'it>n'{t') (10b) 



for the impurity action. The first line in Eq. (jlOp follows 
from Eqs. ()lb|) and (JG]), and the second line results from 
direct insertion of Eqs. © into Eq. px]) . 

Feynman diagrams for pseudoparticle propagators are 
most easily constructed in the grand-canonical ensemble 
with respect to the total pseudoparticle number Q. For 
this purpose it is convenient to switch into the interaction 
representation with respect to a chemical potential term 
XQ. Because Q is a conserved quantity, the value of A 



for i > has no influence on the expectation value of 
physical observables. We choose A to be present only at 
times i < 0, (i.e., on the imaginary part of the contour, 
where t ~ —ir), such that 



am(0 ~ '^m s^PC-^Inr^) (H^) 

alnit) = «]■„ exp(-A Imt), (lib) 

and the grand canonical average can be denoted as 



{m)> 



{e-^^QA{t)) 



Si„ 



(e-^AQ) 



(12) 



Grand-canonical pseudoparticle propagators on the con- 
tour C are defined as 



GLn'it^t') ^ -i{Tca„,it)al,it'))x. 



(13) 



(In the following, propagators are considered to be ma- 
trices in their flavor indices m,m', and the indices will 
be omitted whenever this is not ambiguous.) 

The restricted trace in Eq. ([5]) can be recovered from 
grand-canonical expectation values by means of an ex- 
pansion in powers of the fugacity ( = e"'^'^.— For ob- 
servables which annihilate the Q — state [such as the 
impurity Green's function ([3])], an expansion of Eq. (|12p 
in C. yields 



(Q)a 



{A{t))Q=^ + OiC). 



(14) 



Furthermore, the leading terms of the Green's functions 
(|13p in the fugacity expansion can be obtained in the 
form [cf. Eq. PT|) ] 



g\t,t') = kxit,t')[git,t') + oiO], 



(15) 



where the projected Green's function Q{t,t') is indepen- 
dent of A, and the pre-factor is given by 

fcA(t,t') = e^^""*"""*''[0c(t,i') + 0c(t',t)C]- (16) 

The step function 8c (^,i') is if i is earlier on C than 
t' , and 1 otherwise. In order to obtain a perturbation 
expansion directly in terms of projected quantities, the 
limit A — > CX3 must be taken analytically in all expressions 
below. As a consequence, projected Green's functions 
become the basic objects in the hybridization expansion. 
In addition to the symmetries of the grand-canonical 
propagators?^ projected propagators have a number of 
useful properties, which we list in the following para- 
graph. First, one can show from their definition that 
they satisfy an initial condition 



,(t+,t) 



-iS, 



'mm' 1 



(17) 



when t+ is infinitesimally later on C than t. Furthermore, 
the factor k\ essentially restricts propagation of pseu- 
doparticles to one direction along the contour. In par- 
ticular, the leading order of the product of two Green's 
functions A^ and B^ is given by 

^\t,ii)B^(ti,t') - kx{t,t')A{tM)B{tx,t!) (18) 



if the time arguments t', ti, and t are in cyclic order with 
respect to the arrow in Fig. [U and smaller by a factor C 
otherwise. [In the following, we will use the notation ii < 
t2 -< ... ~< tn to indicate that time arguments ti . . .tn 
are in cyclic order along C, according to the arrow in 
Fig. [T]] Gonsequently, to leading order in C the contour 
convolution of the two functions is given by 



dtA^{t,t)B^{t,t')=kx{t,t') / dtA{t,t)B{t,t') 



C,t'-<t-it 



(19) 
where the integral range on the right hand side must be 
restricted such that t' ~<i -<t. 



B. Pseudoparticle Dyson equation 

Grand-canonical pseudoparticle propagators obey the 
usual Dyson equation with the pseudoparticle self-energy 



/^ = Go + Go *^^*G^ 



(20) 



where [a*b] [t, t') = J^ dt a{t, t) b{t, t') denotes the contour 
convolution, and 



Gi^ra'{t.t') = -^{a„,{t)al,{t')) 



5loo,A 



(21) 



is the bare pseudoparticle propagator (i.e., at zero hy- 
bridization) . The latter satisfies the equation of motion 

[idt-\{t)-E{t)]G^{t,t') = Sc{t,t'), (22) 

where A(i) ~ A on the imaginary part of the contour 
and zero otherwise, and [E(t)]mm' = Smm'Em{t) is a 
diagonal matrix in the flavor indices. We use the no- 
tation of Ref. [l^ for the derivative dt and the contour 
delta- function Sc, i-e-j the latter is defined such that 
J^dtSc{t,t)f{t) = /(t) holds for any function /(t) on 
the C, and dtQc{t,t') = 5c{t,t'). 

Although we will not need an explicit expression for 
the bare projected propagator Go{t,t') in the following, 
it may be a useful illustration to compute it from the 
equation of motion ([2^ and verify that it satisfies all the 
usual symmetries of the contour Green's functions^ and 
the initial condition p7|) . By integrating Eq. (p2|) with 
a periodic or antipcriodic boundary condition for Bose 
and Fermi particles, respectively, we obtain the grand- 
canonical propagator 



6;g\i,t') = -ie^(i«i*-i"i*') 



eyiY>[~i !l, dtE{t)] 



„;3[A+£;(0)] 



x[e^t'+^("%c(t,i') + X0c(t',i)], (23) 

where x = +1 (—1) for Bose (Fermi) particles, and 
Eifi) is the value on the imaginary time axis. Taking 



the limit A — ^ oo in this expressions yields Q^it^t') ~ 
k\{tTt')QQ{t,t'), with the projected propagator 

(24) 

Using Eq. ([22|) . the Dyson equation ean be written in 
differential form 

The corresponding Dyson equation for the projected 
propagators is then derived by inserting Eqs. (|15p . (|16p. 
and (|19p into (|25|) . and taking the limit A — > oo, 



{idt-E{t)\Q{t,i:)- dtJ:{t,i)g{i,t')^0. (26) 

The delta-function on the right hand side has been omit- 
ted, because this equation will be considered only for 

The numerical solution of Eq. (j26p can be performed in 
the same way as the solution of Dyson-like equations for 
real-particle propagators, which is described in detail in 
Ref.[l9l However, the structure of the integral in Eq. ((26)) 
implies an important simplification. In Eq. (|26p . the 
derivative dtg(t,t') is determined entirely by the value 
of g{ti,t') for i' -< ii -< t. For fixed t' , Eq. d^ is thus a 
Volterra integrodifferential equation,— whose numerical 
solution is similar to that of an ordinary differential equa- 
tion with initial condition (J17p . This is particularly inter- 
esting for the initial state, i.e., when all time arguments 
are on the imaginary branch of the contour. Substituting 
the definition (P into Eqs. ^ and ^ yields 



i~dr - A 



E) g^''-' 

- rdf^^'^'{T~f)g 

Jo 



n/3 



A,M / 



<5(r), (27) 



for the grand-canonical version of the Dyson equation, 
and 



{-dr^E)g"{T)- dfl]"(r-f)a"(f)=0, (28) 
Jo 



for the projected Dyson equation. While Eq. ([27)) is a 
boundary value problem and must be solved by Fourier 
transformation [g^'"{P) = ±(?'^'^'(0~) for bosons or 
fermions] , the projected Eq. ()28p is an initial value prob- 
lem [^"(0) = — 1], which is most efficiently solved in the 
imaginary time domain. 



of A can be derived from the standard rules for gen- 
eral quartic interaction terms (see, e.g., Ref . l49i) . Each 
diagram for S^ contains one sequence of pseudopar- 
ticle propagators that connect the two external ver- 
tices (the "backbone"), and possibly additional loops 
of propagator lines, e.g., renormalizations of the hy- 
bridization function (Fig. [5^). To leading order in 
C, the backbone g^{t,tn) ■ ■ ■ g^{t2,ti)g^{ti,t') is given 
hy kx{t,t')g{t,tn)---g{t2,ti)g{ti,t') ffti,...,t„ are in 
cyclic order along C, and smaller by 0{C,) if the vertices 
are not ordered [cf. Eq. ([T5|]. Each closed loop of pseu- 
doparticle propagators contributes an additional expo- 
nentially small factor Q. Thus the diagram rules for the 
projected self-energy T,{t,t') — T,^{t,t')/k\[t,t') ean be 
obtained from the diagram rules for E^ by (i), replac- 
ing pscudoparticlc propagators (|13p by projected propa- 
gators (flS)) . (ii), discarding diagrams with closed loops, 
and (iii), requiring vertices along the backbone to be in 
cyclic order along C. 

For completeness we summarize the final rules for con- 
structing the projected self-energy Yi{t,t'): (i) The nth. 
order contribution to I](t,i') is given by all diagrams 
consisting of 2?i three-leg vertices (Fig. ^p) at times 
to = t',ti,...,f2n = t, of which n correspond to annihi- 
lation operators d (outgoing hybridization line), and n 
correspond to d^ (ingoing hybridization line). The ver- 
tices are labeled according to Fig.[5jD. They are connected 
by one sequence of pseudoparticle lines (solid lines, point- 
ing from t' to i), and n hybridization lines (dotted lines) 
in all possible ways such that the diagram cannot be sep- 
arated into two parts by cutting only one line, (ii) Sum 
over all internal flavor indices, and integrate over the in- 
ternal times ti,...,t2„_i, respecting the cyclic order t' -< ti 
-< ... -< t. (in) Because exactly one fermionic pseudopar- 
ticle operator is attached to each end of a hybridization 
line, the sign of the diagram is (— l)''^-'^, where s is the 
number of crossing of hybridization lines, and / is the 
number hybridization lines that point opposite to the di- 
rection of the backbone, (iv) An overall factor i" msut 
be added. 

The diagrammatic expansion for S ean be resummed 
by replacing bare propagators with interacting propaga- 
tors g, and in turn taking into account only skeleton di- 
agrams, i.e., diagrams in which internal propagator lines 
have no self-energy insertions. Truncation of the skele- 
ton series S]'^'^'''[t/, A] at finite order leads to conserving 
approximations)^ because it ean be derived from the 
Luttinger-Ward functional.— In particular, this fact en- 
sures the conservation of the pseudoparticle number ([7]), 
which is crucial in order to obtain a meaningful approxi- 
mation scheme for nonequilibrium situations. To leading 
order in C, the conservation of {Q)\ implies 



C. Diagram rules for the pseudoparticle self-energy 

Because the local part of the pseudoparticle action 
pop is quadratic, a diagrammatic expansion of pseu- 
doparticle Green's functions and self-energies in terms 



A— >cx) (^ 



J2{~iy^gm{t,t+) (29) 



(a) 



aT^ 



y ■ \ '■., 
i,r?i '■■ — < — ' g •' < ' g • t'vn! 



(b) t,v ■■■<■■■ t,p'=App,{t,t') 

t, m ^<— t, m'= Qmra' {t, t') 

P-><m ={m{t)\dp\m'{t)) 



(c) 






iii) 




< 




V) 
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,.-••■••■ 


""?^\ 


■\ 



< • g ■ < 



■ '' < ' •' < ■ •' < : '' <'• c 



vi) ^.^■■■■'p'--'y^-':~" ■■■■.. 
•^—« — ■ •' « '"' < — ' ■ g '• ■ g 



(d) 




FIG. 2: (a) A 5th-order contribution to the self-energy 
^mm' (*i*')i consisting of Green's functions 5o (solid lines) 
and hybridization functions A (dotted lines), (b) Building 
blocks of the diagrams in the hybridization expansion. The 
pseudoparticle line (solid line) corresponds to the interacting 
propagator Q in skeleton expansions, and to Qo otherwise, (c) 
AU diagrams for E'""''[^, A] up to third order. In the topolo- 
gies indicated, the hybridization line can point in any direc- 
tion, which gives 2, 4, and 4x8 diagrams in first, second, 
and third order, respectively, (d) Factorization of third order 
diagrams [(iii)-(vi) in (c)] by separating out the vertex part. 



where t+ is infinitcsimally later on C than t, and (—1)'" = 
±1 if m corresponds to Bose or Fermi particles, respec- 
tively. These relations provide a good check for the nu- 
merical implementation. 

All skeleton diagrams up to third order are displayed 
in Fig. [2t. The self-consistent strong-coupling expansion 
has been proposed long ago^"— as an approximate so- 
lution for the Anderson impurity model. Kuramoto^ 
coined the term non-crossing approximation for the low- 
est order, i.e., keeping only the first diagram in Fig. [5}:. In 
the present work we use the skeleton series up to third or- 
der as an impurity solver within nonequilibrium DMFT, 
and compare the results to CTQMC (Sec HV)) . 



D. Diagram rules for the impurity Green's function 

In general, expressions for observables in the impu- 
rity model can be derived from the grand potential 



fix ~ — /5~^ logTr[C'^Tce'^''°p]. Because diagrams for the 
correction AU\ = fix — n\{A = 0) contain at least one 
closed loop of pseudoparticle lines, Ail\ is proportional 
to C for A — > oo. The leading order in (,, 



Q 



-1 



lim _logTr[C^Tce 



o5i„ 



(31) 



is obtained by adding to the local contribution 51(A = 0) 
all diagrams of fl\ which contain only one loop, in which 
Q^ is replaced by G, and where integrals over the internal 
vertices are restricted such that the vertices are in cyclic 
order. (See the analogous argument for E in Sec. lIIICl ) 
Using Eq. (|14l) . the impurity Green's function ([3]) 
is given by G{t,t') = limA->oo G'^(i,i')/('3)A, where 
G^p, {t, t') = -i{Tcdp{t)dl, (i'))A- It can thus be obtained 
from the derivative [cf. Eqs. dH), dTU]), dMl), and ^] 
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Lrpp' [Z, T J — -. 



Sfl 



Q5kp>p{t',ty 



(32) 



Diagrams for G{t,t') (in terms of the projected pseu- 
doparticle Green's functions) are therefore constructed 
by removing one hybridization line from the diagrams 
for n (Fig. [3]) . Note that a diagram for H generally has 
a symmetry factor 1/5* 7^ 1, where S is the number of 
topologically equivalent ways to label the vertices. The 
symmetry factor disappears in the expansion of G, be- 
cause for a diagram with symmetry factor 1/S there are 
S ways to remove a hybridization line which lead to the 
same diagram for G. The series for G can thus be re- 
summed in the same way as the series for E, i.e., by 
keeping only skeleton diagrams for G, and replacing Go 
with G- Equation ([5^ then holds also for the skeleton 
expansioir. 



Gllf[G,A]{t,t') 



psn'''^'[G,A] 

Q SAp'p{t',t) '' 



(33) 



where r2'^'^'''[G, A] is the Luttinger Ward functional, i.e., 
the skeleton expansion for fl in terms of the fully inter- 
acting (projected) propagators G- To design an approx- 
imation for G^'^'^^[G,A] which is consistent with a given 
approximation of E one must truncate both 0'^'^'''[C7, A] 
and E'^'^'''[t/, A] at the same order. 

The final rules for G^'^'^^[G, A] read: The nth order con- 
tribution consists of a loop of projected pseudoparticle 
propagator lines (Fig. [2}d) which connects 2n vertices (n 
annihilation operators, n creation operators). One d- 
vertex (time t, A-line labeled p) and one d^-vertex (time 
t', A-line labeled p') are external vertices. The inter- 
nal vertices are connected by hybridization lines such 
that no internal line has a self-energy insertion. Sum 
over all internal flavor indices, and integrate over inter- 
nal (contour) time variables respecting the cyclic order 
of ii . . .t2n along the contour. Add a pre-factor z". To 
determine the sign of a diagram D, reinsert the A-line 
between the external vertices. This recovers the diagram 
D' in the expansion of f2 from which the diagram D is 




FIG. 3: All diagrams for G^'^^'ia.A] in first order [diagram 
(i)], second order [diagram (ii)], and third order [diagram (iii)- 
(vi)]. Internal hybridization lines can point in any direction, 
which gives f , 2, and 4x4 diagrams in first, second, and third 
order, respectively. 



derived. The sign kjj of D is given by the sign of D', 
i.e., kd = {~^y~^^ : where s is the number of crossings 
of hybridization lines, and / is the number of hybridiza- 
tion lines that go in opposite direction to the string of 
pseudoparticle lines that is obtained if the loop is cut 
at an arbitrary fermionic propagator line. All skeleton 
diagrams for G up to third order are shown in Fig. [3] 



The numerical effort of the evaluation of diagrams is 
mainly determined by the contour integrals over the in- 
ternal vertices. Using Monte Carlo for the evaluation of 
the integrals for higher order diagrams will suffer from 
a sign problem. We use a quadrature formula for the 
integrals which is based on equidistant discretization of 
the contour C. The nth order diagrams have 2n — 2 in- 
ternal integrals (cf. Figs. [2t and [3]), which have to be 
evaluated for each combination of the two external time 
variables. (Nonequilibrium correlation functions depend 
on both time arguments separately). This seems to im- 
ply that the numerical effort for the evaluation of E and 
G scales with the number N of mesh points like N"^ , N^, 
and N^ for first, second and third order, respectively. 
(The effort for the solution of the Dyson equation, which 
is essentially a matrix inversion on C, scales like N^.) 
However, one can reduce the effort for the evaluation of 
the third-order diagrams to N^ by factorizing out a ver- 
tex part F with two internal integrals and three external 
variables (Fig.[2Ji); E and G can then be computed from 
F with only two additional integrals. Since F does not 
have to be stored in memory, this way of evaluation is 
more efficient than performing four internal integrals. 



E. Numerical implementation 



IV. COMPARISON TO CTQMC 



Before presenting first benchmark results for the self- 
consistent hybridization expansion, we would like to 
make some remarks on the numerical implementation. 
First of all, we note that the number of possible labelings 
for the internal flavor indices is usually quite restricted. 
As an example, consider the single-impurity Anderson 
model with the four basis states |0), jcr) = rf^|0), and 

in) = 44io)' ^1- = ud\d^did^ ~ fiid\d^ + djd^), 

and Shyb = -ijrdtdt' J2a diit)A^{t,t')d^{t'). The ma- 
trix elements (Pq) are nonzero only for the combinations 
F^Q = 1 and F^. ^ = a. Furthermore, Green's func- 
tions are diagonal in pseudoparticle flavor, because both 
the interaction part and the hybridization function are 
diagonal in the occupation number basis. The second- 
order digram for Emm(t,t'), e.g., then only allows eight 
possible labelings for the three internal Green's function 
lines, which are (t, 0, 1) and (|, 0, t) for m =t|, (0, a, ti) 
and (t4-, 0-, 0) for m = a, and (t, tl, i) and (I, tl, t) foi' 
m = 0. 

To obtain a self-consistent solution, the hybridization 
expansion is evaluated by iteratively solving the Dyson 
equation (pS)) for Q{t,t'), and evaluating the integral ex- 
pressions for E(t, i'). However, the real-time version of 
the expansion can easily be implemented in a slightly 
simpler way that exploits the causal structure of the 
equations. If we have computed G{t,t') for Re(i), Re(i') 
< imax, then ^(iniax + Ai, t') and g{t, i^ax + A) can be 
obtained from the above mentioned iteration in only one 
or two steps by starting from a polynomial extrapolation 
of g{t, t'). This amounts to a stepwise propagation of the 
solution on the imaginary branch of C to real times. 



A. Interaction quench in the Hubbard model 

Nonequilibrium DMFT for the interaction quench in 
the Hubbard model provides a perfect framework to 
benchmark the perturbative impurity solver. The Hub- 
bard Hamiltonian 



m) = E ^'j-^s- + ^(^) E ("''t-s) ( 



nii' 



(34) 



describes fermions of spin one half which hop on a lattice 
with hopping amplitude Vij and interact with a repulsion 
energy U on each site. To perform an interaction quench, 
the system is prepared in a thermal equilibrium state at 
temperature T = 1//3 and interaction U{t < 0) = C/q for 
times t < 0, and the interaction is suddenly switched to 
a new value U{t > 0) = U at t = 0. 

The DMFT equations for the interaction quench have 
been explained in detail in Refs. [l^ andfl^. In the follow- 
ing we assume that the hopping matrix Vij has a semiel- 
liptic local density of states 



p{e) = ^4-{e/Vy/2TTV 



(35) 



(with half-bandwith 2V), and we focus on the paramag- 
netic state at half-filling. DMFT then reduces to a set 
of two self-consistent equations >^'i^ (i) The local lattice 
Green's function ^ must be determined from the single- 
site action ([T]) , where the index p now labels spin a =t, 4-, 
and the local Hamiltonian is given by 



i7loc(t)-f/(t)EK- 5)^-5) 



(36) 
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FIG. 4: Double occupancy deq{P,U) in the thermal equilib- 
rium state. The impurity solver is either CTQMC, or the self- 
consistent hybridization expansion up to first (NCA), second 
(OCA), and third order, (a) /3 = 5. (b) /3 = 10. (c) /? close to 
the critical temperature Tc of the first-order metal-insulator 
transition. 



(ii) The hybridization function is determined by the self- 
consistencj*^ 



K,{t,t')^V^G,{t,t'). 



(37) 



The hopping 1/ = 1 is used as an energy unit, and times 
are measured in units of the inverse hopping [h = 1). 

Below we solve these DMFT equations by means of 
the self-consistent hybridization expansion and compare 
to results from CTQMCJ^ii^ In particular, we focus on 
the time-evolution and the thermal equilibrium value of 
the double occupancy per site, d{t) = {ni^riii), which is 
a local observable and can thus be measured directly in 
the impurity model, i.e.. 



dit) = iQ-'gn{t,t+), 



(38) 



where Q-^i is the propagator for doubly occupied sites, 
and f^ is infinitesimally later on C than t. 



B. The initial equilibrium state 

Figure ID shows the double-occupancy deq{l3, U) in the 
thermal equilibrium state at interaction U and inverse 
temperature /3. At large enough temperature, deq{T, U) 
decreases smoothly as a function of U (Fig. |3^). As T 



is lowered, the curves bend strongly around U = 4.5 
(Fig. Hb), indicating a narrow crossover between metal- 
lic and insulating behavior. Below a critical temperature 
Tc, the transition between metal and insulator becomes a 
first-order phase transition (Fig.|4]:). For the semielliptic 
density of states, the endpoint of this Mott transition is 
located at Ur = 4.7 and T. = 0.055^ 



In agreement with recent results based on a Monte 
Carlo sampling around NCA^ we find that NCA can 
reproduce the CTQMC results only deep in the insulat- 
ing phase. However, already the lowest order correction 
to NCA, i.e., OCA, very well accounts for the nonlinear 
behavior of deq{U,T) in the crossover regime, and the 
third order in the self-consistent hybridization expansion 
almost quantitatively recovers the CTQMC results even 
close to the critical point (Fig. Ht). The location of the 
critical endpoint in the phase diagram is in good agree- 
ment with previous QMC results)^ If we estimate Tc 
from the smallest /3 for which we can detect hysteretic 
behavior in deq{/3,U) (this gives actually a lower bound 
for Tc), we find Tc > 1/19 « 0.052 for the 3rd order 
and Tc > 1/26 « 0.038 for OCA (Fig.H;). NCA, on the 
other hand, does not display singular behavior in this pa- 
rameter regime. As usual, the convergence of the DMFT 
equations slows down close to the critical point, and it is 
thus hard to get precise numbers for Tc and Uc- 

The order-by-order convergence of the self-consistent 
hybridization expansion is also evident from the local 
Green's function G{t), both in the crossover regime 
(Fig. [5ji) and in the insulating phase (Fig. \5jp). From 
the value G{(3/2) one can see that NCA overestimates 
the insulating nature of the solution. This fact reflects 
a well known deficiency of NCA: The Kondo tempera- 
ture Tk for the Anderson model comes out correct within 
NCA for [/ = oo, but it is severely underestimated for fi- 
nite interaction U. This problem can be cured by taking 
into account certain vertex corrections, which correspond 
to summing up higher order terms in the self-consistent 
expansioui^i^S Our results show that the third order is 
sufficiently accurate in a wide parameter range covering 
the insulating phase and the crossover regime, even close 
to the critical point. 

Another known deficiency of the NCA is that the 
Fermi-liquid in the Anderson impurity model for T <^ Tk 
is not correctly describedi^ Because this problem cannot 
be cured by taking into account finite order diagrams 
in the hybridization expansiouj^S one would expect that 
even the third order will yield wrong results in the metal- 
lic phase at low enough temperature. Empirically, we 
find a slow-down of the convergence as the temperature 
is lowered in the metallic phase, and we have not sys- 
tematically studied the breakdown of the truncated self- 
consistent hybridization expansion deep in the metallic 
phase. 
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FIG. 5: Matsubara component ((5} of the local Green's func- 
tion Gcr{T) [Eq. [3], which is independent of a in the param- 
agnetic phase. The impurity solver is either CTQMC, or the 
self-consistent hybridization expansion up to first (NCA), sec- 
ond (OCA), and third order, (a) U = 4.5, /3 = 20 (crossover 
regime), (b) U = 6, P = 20 (insulating phase). Note that the 
CTQMC results in the insulating phase are only accurate for 
values larger than the "noise floor" of about 10"'^. Panels (c) 
and (d) show the same data as (a) and (b), respectively, but 
plotted for small r on a linear scale. 
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FIG. 6: Double occupancy d{t) after an interaction quench 
from f/o to U. (a) and (b): initial states in the crossover 
regime [Uo = 3, /? = 5], (c) and (d): initial states in the 
insulating phase [Uo = 5, /3 = 5]. 



MODULATION SPECTROSCOPY ON THE 
MOTT INSULATING PHASE 

A. Introduction 



C. Time evolution of the double occupancy 

To test the accuracy of the strong-coupling expan- 
sion for noncquilibrium problems we compute the time- 
evolution of the double occupancy after an interac- 
tion quench. Due to the dynamical sign problem, 
weak-coupling CTQMC calculations for interacting ini- 
tial states^ can be performed only to relatively short 
times tmax! which mainly depend on the final interac- 
tion U{t > 0). However, for those imax which are ac- 
cessible with CTQMC, the comparison with the strong- 
coupling expansion reveals a similar trend as for ther- 
mal equilibrium states (Fig. [6|) : For quenches from the 
crossover region to larger interaction both the initial state 
and the time evolution is not correctly described within 
NCA, whereas OCA is more reliable, and the third order 
calculation recovers the CTQMC results almost quanti- 
tatively (Figs. |B^ and b). While Figs. [B^ and b show 
the longest times accessible with CTQMC, the OCA and 
the 3rd order calculations can be carried to substantially 
longer times. For a quench from the insulating phase 
to smaller interaction, NCA is better suited to describe 
the initial state, but differences to CTQMC become more 
pronounced during the time evolution (Figs.[BJ: and d). 



To illustrate the capability of the approach, we are now 
going to present an application of the strong-coupling 
hybridization expansion in a parameter regime where 
the weak-coupling CTQMC approach would be numer- 
ically too expensive. Our aim is to compute the re- 
sponse of the double occupancy in the Mott insulator 
to a time-dependent change of the interaction U or the 
hopping amplitude V. Such an experiment, with a peri- 
odic change of the hopping, was originally proposed by 
KoUath et ali^ as a new type of spectroscopy for ultra- 
cold gases in optical lattices without direct analogon in 
solid state physics. In the meantime, the technique has 
been used as an experimental probe for the detection of 
the Mott insulating phase of ultracold ^°K-atoms in an 
optical latticei^ 

In the experiments, the hopping amplitude V{t) = 
Vb[l + acoii{u)t)] is modulated sinusodially over several 
tens of periods 2ti / uj ^^■ ' '^'^ Apart from an oscillating com- 
ponent dosc{t) with zero time-average, the double occu- 
pancy d{t) rises linearly in time for small times, and sat- 
urates within a timescale Tgat- In the Mott insulating 
phase, the modulation spectrum, i.e., the magnitude of 
the response as a function of the frequency, has a peak 
when uj is approximately at resonance with the energy U 
that is needed for the creation of a doublon-holon pair. 
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and a gap at w = 03i^ 

Because the modulation strength can be quite large, 
many aspects of those experiments can only be under- 
stood by means of a nonequilibrium formalism. This 



0.16 



is certainly the case for the saturation time r,, 



55 



and 



for the nonequilibrium quasisteady (time-periodic) state 
which the system is in once it has saturated. Even when 
averaged over time, such a state might have properties 
which do not resemble any thermal equilibrium state of 
the system. Nonequilibrium DMFT can be used to re- 
solve some of these issues. In the following, we demon- 
strate that DMFT yields a modulation spectrum which is 
in agreement with a recent investigation based on slave- 
boson mean-field theory;^ and similar to what has been 
obtained in timc-dcpcndcnt density-matrix renormaliza- 
tion group calculations for the one-dimensional Hubbard 
modeli^ Furthermore, we will show that a slightly differ- 
ent modulation procedure, namely a quench of the hop- 
ping or the interaction by a few per cent, provides another 
probe which is sensitive to the Mott transition and the 
metal-insulator crossover at higher temperatures. 



B. Periodic modulation of U 

In the following we consider the Hubbard model ([M)) 
with a time-dependent interaction 



U{t > 0) = Uo[l + acos{ujt)], 



(39) 



where a is the relative modulation strength (a < 1). For 
times t < 0, the system is prepared in an equilibrium 
state at interaction U{t < 0) = Uq and temperature T = 
1//3. The model is treated within nonequilibrium DMFT, 
assuming a semielliptic density of states (PSJ) , such that 
the self-consistency is given by Eq. (|37)) . The energy scale 
is fixed by the quarter bandwidth V = I. Because we will 
restrict the investigation to insulating states and to the 
crossover regime, we will mainly use OCA as an impurity 
solver. 

Experimentally, the hopping is more easily tunable 
than the interaction, because the former is strongly af- 
fected by a change of the optical lattice depth. Our 
description ([55]) is nevertheless justified, because mod- 
ulation of U and V are equivalent from a theoreti- 
cal point of view, and only the ratio U/V matters. 
The Hubbard model H{t) with time-dependent inter- 
action (j39p and time-independent hopping Vo can be 
mapped onto an equivalent model with time-independent 
interaction C/g a-nd periodic hopping V{t') = Vo/{l 4- 
acos[u}t{t')]}, where i(i') is the inverse of the transfor- 
mation t'{t) = t + a sm{ojt) / uj . To lowest order in a, 
modulation of U{t) ~ Uo[l + acos{ujt)] and V{t') ~ 
Vb[l — acos(a;i')] are thus equivalent, although higher 
harmonic terms oc a" cos{nu!t') appear in V{t') for large 
a. The equivalence can be established by a simple 
change of time-variables in the time-evolution opera- 
tor U{t,0) = Ttexp[-iJl^diH{i)] from t to t', which 
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FIG. 7: (a) Double occupancy d{t) during a periodic modu- 
lation of the interaction around Uo = 7 {13 — 10, 5U — alio — 
2). The frequency is below the resonance peak {u — 1.5), close 
to the resonance {u — 7), and above the resonance (oj — 10). 
Bold solid lines indicate the average dav{t) over one period 
[t — n/bj, t + n/ijj]- OCA was used as an impurity solver, (b) 
Modulation spectrum for Uo — 7 and /3 — 10, using either 
NCA or OCA as an impurity solver: For large amplitude 
SU = 2, the data have been obtained from dav{to) at given 
time to = 10 (open symbols); for small amplitudes 5U = 0.5, 
the slope {d/dt)dav{t) in the interval 8 < i < 10 is plotted 
(full symbols). The curves have been scaled with a constant 
factor to match their peak heights, (c) Averaged double oc- 
cupancy davit) at Uo = 7 for various modulation amplitudes 
SU = aUo, plotted against tSU'^ . The three curves labeled 
NCA almost fall on top of each other. 



yields U'{t',0) = U{t{t'),0), where U'{t',0) is the time- 
evolution operator for the model H'(t') with periodically 
modulated hopping V{t'). 

The typical time-evolution of the double occupancy 
d{t) during a periodic modulation of U is displayed in 
Fig. [7^. The effect of the modulation vanishes for w — 5- 0, 
and is largest close to tu ~ U. The anharmonic behavior 
in d{t) for small frequencies w is due to a rather large am- 
plitude SU = alio. After removing the oscillating compo- 
nent dosc{t) by taking an average dav{t) = Jf_^/2 dt' d{t') 
over one period r = 27r/a;, one can clearly identify an ini- 
tial linear increase with a slope r(w), followed by a trend 
towards saturation at longer times. 

The rate r(a;) can be obtained from second order 
time-dependent perturbation theory, i.e., F oc a^ for 
a ^ 1^2^ (Linear response yields only an oscillating 
contribution.) Since the resulting Fermi's golden rule ex- 
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pression for F depends only on the equilibrium state of 
the system at U = Uor^ it would be best to measure 
r(a;) directly in the limit a -t> 0. To extract the slope 
r(a;), one has to go to small amplitudes a anyway, be- 
cause only then does the linear region extend over many 
oscillation periods. On the other hand, experiments are 
performed at quite large a in order to obtain a good sig- 
nal to noise ratio, and the magnitude of the response is 
defined by the value of the double occupancy dav (t) at 
given time Iq. Usually, ip is chosen so large that dav{t) 
is no longer linearly increasing at i = io- 

In Fig. [TId we compare two ways to quantify the re- 
sponse: Either (i), the increase of the double occupancy 
davito)—d{0) is measured at a given time to for large mod- 
ulation amplitudes as a function of frequency, or (ii), the 
slope F(w) is obtained from a linear fit to the initial in- 
crease for smaller amplitudes a. Both approaches give a 
modulation spectrum with a peak at oj « C/, and a gap at 
w = 0. Similar to the findings of the previous section, the 
NCA solution slightly overestimates the insulating nature 
of the solution compared to the more reliable OCA. Like 
in the one-dimensional casc;^ our data show that the lo- 
cation of the peak at w « t/ is not considerably shifted if 
the measurement is no longer performed at a <C 1. The 
gap, on the other hand, is most reliably extracted from 
the second approach (ii). In the following we will only 
focus on the peak, and not investigate the low frequency 
weight in detail. 

Figure [5^ displays the peak in the modulation spec- 
trum for various values of C/q . For the semielliptic density 
of states, the first-order phase transition line terminates 
at Uc ~ 4.7 and Tc « 0.055^^ and the zero-temperature 
transition is located at Uc2 ~ 6. Hence the data in 
Fig. [8l which are computed at T ~ 0.1, correspond to 
a cut through the crossover region of the metal-insulator 
phase diagram. The peak in the modulation spectrum 
is clearly visible all throughout the insulating phase and 
the crossover regime between metal and insulator. Its po- 
sition a;inax(C^) scales linearly with Uq in the insulating 
phase (Fig.[5}D), while in the crossover regime we find only 
a weak dependence on Uq. This finding is in good agree- 
ment with results for r(a;) from slave-boson mean-field 
theory?^ where peaks around uj ^ U and cu sa Uc2 are 
predicted for the insulating phase and the metallic phase, 
respectively. In the slave-boson approach, these spec- 
tral features arise from excitations between the Hubbard 
bands in the insulator, and between pre-formed Hubbard 
bands in the metallic phase. 

Since our calculation is performed at temperatures 
above Tc and the pre-formed Hubbard bands shift with 
U, the almost kink-like Wniax(C^) seems very remarkable. 
It should be investigated whether the third order in the 
expansion leads to a smoothening of this feature. This 
will require quite some numerical effort (although it is 
still feasible using small-scale parallelization) , and it is 
thus left to future work, which should involve a more 
realistic setup. 

Another interesting topic is the saturation of dav (t) at 
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FIG. 8: (a) Modulation spectrum at /3 = 10 and various 
interactions Uq, obtained from tlie increase dav (to) — d{0) at 
to — 10 {SU = 2). (b) Location of the maximum ojmax of the 
spectra in (a), plotted against Uq. 



large timesi^ A scaling of the time-axis with a^ indicates 
that the saturation time Tsat behaves like Tgat ct a~^ in 
the insulating phase for frequencies u) — U (Fig. [Tj;), 
while the saturation value dav{t — > cxd) does not depend 
sensitively on 6U. However, due to the long times needed 
to reach saturation at small SU, this result has so far only 
been computed using NCA as an impurity solver, which 
is reliable only deep in the insulating phase. A depen- 
dence Tsat oc a^'^ would be consistent with an incoherent 
pumping mechanism into a doublon-holon continuumj^ 



C. Quench-like modulation of U 

Modulation spectroscopy for the double occupancy is 
in principle not restricted to the periodic modulation 
Eq. (j39p . In particular, a small interaction quench from 
U{t < 0) = t/o to U{t > 0) = Ua + SU, or an equiva- 
lent quench of the hopping, can be viewed as a modula- 
tion experiment as well. In order to extract a frequency- 
dependent response signal, we define the Fourier trans- 
form 

d{oj) = Re / "Ii e^"* [d{t) - d(t„,ax)] , (40) 

Jo 

where imax is the maximum time reached in the simula- 
tion, and d{t) = {l/Z)TT[e-^"°eJ"*de-'"^] is the time- 
dependent double occupancy. Using first-order perturba- 
tion theory one obtains 



d{oj) = w-i Imx(w) + 0[{6U/3f 



(41) 



where x(w) = J^ dseA'^'^^^'^^xi^) is the Fourier trans- 
form of linear response function 



Xis) = -^-Tr(e-^^[e'^^■^(ie-^^^d]) (42) 
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FIG. 9: (a) The double occupancy after quenches from U — 
Uo to U = Uo - SU, with 517 = 0.5, /3 = 10, and Uo = 
7, 6, 5.6, 5, 4.5. OCA has been used as an impurity solver, (b) 
Fourier transform (|40|) of the date in (a). 



of the final Hamiltonian H{t > 0). Equation (PT|) holds 
for uj ^ and tmax ^^ oo, while evaluation of d{uj) at 
imax < corresponds to a broadening of xi^) on the 
scale 1/tmax- The time-independent expansion param- 
eter PSU in Eq. (|4T|) is obtained by performing a per- 
turbation expansion of the initial state density matrix 
exp[—l3H{t < 0)t] instead of the time-evolution operator 
exp[-iH{t> 0)t]. 

In principle, xi'-^) could be computed directly from 
equilibrium DMFT. However, the equilibrium imaginary 
time formalism requires an analytical continuation, while 
the nonequilibrium calculation gives direct access to fre- 
quency or time-dependent quantities. Furthermore, in 
experiment the values of U would have to be changed by 
at least a few per cent, such that it is unclear whether 
Eq. (PT|) is still appropriate to describe the response. 

After an interaction quench from the noninteracting 
ground state to the insulating phase in the Hubbard 
model, the double occupancy d{t) relaxes through a se- 
ries of oscillations at approximate frequency U^^ which 
correspond to the well-known collapse and revival os- 
cillations in the limit U 3> V^ Similar oscillations be- 
come visible after small quenches within the insulating 
phase ([/+ = 6.5,5.5,5 in Fig. [S]); their Fourier trans- 
form leads to a broad peak around uj k, U (Fig. [91d). 
The series of quenches displayed in Fig. [S] corresponds to 
a scan through the phase diagram at constant tempera- 
ture T = 0.1 above the critical temperature of the Mott 
transition. A crossover between metallic and insulating 
behavior is thus expected between C/ « 4 and C/ « 5. In 
agreement with this, the gap in diuj) at w = disappears 
between U = 5 and U = 4.75. Although the absolute 
changes of the double occupancy are small, these results 
suggest that the double occupancy after an interaction 
quench may be used as probe of the metal-to-insulator 
crossover in the Hubbard model. 



VI. CONCLUSION 

In this paper we have presented a self-consistent di- 
agrammatic strong-coupling expansion on the Keldysh 
contour, which can be used to solve the Anderson impu- 
rity model in rather general nonequilibrium situations. 
The main purpose of this study was the development of 
an impurity solver for nonequilibrium DMFT which can 
cover the regime of large interactions and relatively long 
times. By comparing results from our strong-coupling ex- 
pansion to numerically exact weak-coupling CTQMC for 
the Hubbard model, we found that the strong-coupling 
expansion is a good candidate to fulfill these require- 
ments: While it fails in the metallic phase and at very low 
temperatures, and while the first order or non-crossing 
approximation is correct only deep in the insulating 
phase, an important correction arises already in the sec- 
ond order (OCA). In the insulating phase and in the 
crossover regime the latter gives quite reliable results, 
which can be brought into almost quantitative agreement 
with CTQMC by going to the third order of the expan- 
sion. 

Although the numerical effort for the evaluation of 
the diagrams rises considerably with the expansion or- 
der, even calculations up to third order can be carried 
to substantially longer times than CTQMC, because the 
latter suffers from an exponential increase of the com- 
putational cost due to the dynamical sign problem. We 
thus believe that the strong-coupling expansion will al- 
low to extend nonequilibrium DMFT investigations into 
the parameter regime of rather strong interactions and 
not too low temperatures which was so far not accessi- 
ble within weak-coupling CTQMC. It is precisely this 
parameter regime which is relevant for many experi- 
ments with cold atomic gase o^^'^^ and pump-probe spec- 
troscopy. The broad range of possible applications of 
DMFT in this field includes the excitation of the Mott 
insulating phase by a short laser pulse (similar to the 
experiments in Refs. y and|j), the respose of the Mott 
insulator to very strong electrical fields that might lead 
to a dielectric breakdown;^ or an extended investigation 
of the interaction quench in the Hubbard model J^ 

In the last part of the paper we have used the self- 
consistent hybridization expansion as an impurity solver 
within nonequilibrium DMFT in order to study the gen- 
eration of doubly occupied sites in a Mott insulator by 
a time-dependent variation of the interaction or hopping 
strength. In agreement with previous investigations r^*^ 
the modulation spectrum was found to have a gap at 
a; = and a pronounced maximum aX u Ki U . The max- 
imum persists in the crossover regime, but its location is 
then no longer proportional to U . Furthermore, we have 
studied the double occupancy d{t) as a function of time 
after small interaction quenches. In the crossover regime, 
the behavior of d{t) is drastically changing, which is most 
clearly evidenced by the disappearence of the gap in the 
Fourier transform of d{t) . Although the absolute changes 
of the double occupancy are small, its time-evolution 
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may thus yield a sensitive probe of the metal-to-insulator 
crossover in the Hubbard model. 

An interesting next step would now be to repeat simi- 
lar calculations in a more realistic setup, e.g., on a cubic 
lattice. A careful comparison of results from the self- 
consistent hybridiation expansion up to second and third 
order, and from CTQMC (for small times) will allow to 
make definite experimental predictions, such as for the 
modulation spectrum in the insulating phase and the 
crossover regime, and for the saturation behavior of the 



double occupancy. 
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